AEROSPACE  RESEARCH  •  AER00YN4MICS  •  PROPULSION  •  STRUCTURAL  OYHAMICS  •  ELECTRONIC  SYSTEMS  ANO  INSTRUMENTS  •  COMPUTER 


2 


A 


> 


* \ RESEARCH 
ENGINEERING 
PRODUCTION 

l  •* 

.  .* 


TECHNICAL  REPORT  NO.  549 


COMPUTER  PROGRAMS  FOR  A  REACTIVE 


TURBULENT  BOUNDARY  LAYER  - 


AIR  VERSION 


Bv  B.  Bellow 


CLEARINGHOUSE 

FOB  FEDERAL  SCIENTIFIC  AND. 
technical  information 


Microfiche 

^-7r 


Ml 


September  15,  1965 


GENERAL  APPLIED  SCIENCE  LABORATORIES,  INC. 

MERRIQK  and  STEWART  AVENUES.  WEST^URY,  1_.  I:,  N.V.  ■  (516)  ED  3-69p0 


BEST 

AVAILABLE  COPY 


Project  8030/1220 


Total  No.  of  Pages  -  v  and  81 
Copy  No.  (  9  )  of  135 


TECHNICAL  REPORT  NO.  549 
COMPUTER  PROGRAMS  FOR  A  REACTIVE 
TURBULENT  BOUND  .RY  LAYER  -  AIR  VERSION* 

By  B.  Bellow 


Prepared  for 

Advanced  Research  Projects  Ac  icy 
Washington  25,  D.  C. 


Under  Contract  SD-149,  Supplement  #5 
ARPA  Order  No.  395 
Project  Code  3790 
"Ballistic  Reentry  Studies" 


Project  Engineer  -  W.  Daskin 
Code  516  -  ED  3-6960 


Prepared  by 

General  Applied  Science  Laboratories,  Inc. 
Merrick  and  Stewart  Avenues 
Westbury,  L.  I.,  New  York 


September  15,  1965 

k 

This  research  was  sponsored  by  the 
Advanced  Research  Projects  Agency 


Approved  by: 


Antonio  Ferri 
President 


11 


SUMMARY 


Computer  programs  for  the  calculation  of  properties  within 
a  reactive  compressible  turbulent  air  boundary  layer  on  a  flat 
plate  are  described.  The  pressure  is  constant  throughout  the 
boundary  layer.  The  partial  differential  equations  for  energy 
and  species  mass  conservation  are  solved  with  arbitrary  initial 
conditions  by  a  finite  difference  technique.  A  variable  wall 
temperature  boundary  condition  may  be  used.  The  boundary  con¬ 
ditions  at  the  edge  of  the  boundary  layer  are  constant  with 
respect  to  the  axial  coordinate.  The  partial  differential  equa¬ 
tions,  which  describe  a  two-dimensional  diffusing  flow  may  be 
coupled  to  a  system  of  equations  describing  a  one-dimensional 
finite-rate  air  chemically  reacting  flow,  and  hence  may  be  used 
in  the  numerical  treatment  of  the  two-dimensional  r eacting 
boundary  layer. 


iii 
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TECHNICAL  REPORT  NO.  549 


COMPUTER  PROGRAMS  FOR  A  REACTIVE 
TURBULENT  BOUNDARY  LAYER  -  AIR  VERSION 


By  B.  Bellow 


INTRODUCTION 

This  report  describes  an  IBM-7094  computer  program  written 
in  the  FORTRAN  II  language  to  calculate  properties  within  the 
turbulent  boundary  layer,  with  air  chemistry.  The  analysis  is 
described  in  Ref.  1.  There  are  three  versions  of  this  pro¬ 
gram,  reflecting  the  substructure,  reference,  and  sublayer 
hypotheses. 

Section  I  of  this  report  will  describe  the  details  common 
to  all  versions.  Sections  II  and  III  will  outline  those 
features  peculiar  to  the  substructure,  reference,  and  sublayer 
versions.  The  major  differences  in  the  three  versions  are  the 
computation  of  the  parameter  d/d\(^n  <J)  and  the  input  format 


for  execution  on  the  IBM  7094. 
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I.  TURBULENT  BOUNDARY  LAYER-AIR  CHEMISTRY 
A.  Basic  Equations  Used 

The  program  solves  two  partial  differential  equations 
for  the  dependent  variables,  stagnation  enthalpy  ratio,  "G", 
and  species  mass  fractions,  "Y^." »  as  functions  of  the  two 
independent  variables  X  and  0.  These  equations  are: 


3G  _ 

a 

2 

~  U 

u  oG  e 

[.  1 

~  ilHll 

U 

P  *  S0  2h 

e  e 

~J 

e 

*  dx  (in  a) 


dG 

bij) 


(1) 


_d_ 

&Y.  1 
u  k 

,  d  . . 

_ 

30 

Icn 

a0 

(2) 


[  k  =  1  to  7  in  Eq. (2) 
referring  to  species 
0,N,N0,C>2  ,02  ,Na  ,N0+] 

Explanation  of  symbols: 

Constants  in  coefficients 

P  -  Prandtl  number 
e 

S  -  Schmidt  number 
e 

u^  -  Reference  velocity 
h^  -  Reference  enthalpy 
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Parameters  that  are  dependent  on  \  and 


u  =  u/u  ,  G  =  - — 
e  h 

e 


(3) 


«  =  =“  {(7  +1  V  aj  < ln  [9(C.x>  ■>} 


(4) 


s 

Equations  (1)  and  (2)  are  converted  to  difference 
form  and  solved  by  a  tri-diagonal  matrix  method  in  subroutine 
STEP.  (See  Section  I-B  for  method  of  solution.)  The  species 
equation  [Eq.  (2)]  is  solved  first,  followed  by  the  energy 
equation  [Eq.  (1)].  A  two-dimensional  grid  of  lattice  points 
in  the  (£,$)  plane  is  constructed  as  shown  in  Fig.  1. 

The  properties  of  these  mesh  are  identified  by  suit 
able  subscripts  and  superscript  as  described  below: 

Subscript,  i  -  mesh  points  in  ip  direction 
k  -  scans  over  7  species 
Superscript,  n  -  mesh  points  in  £  direction. 


EXPANSION 

REGION 


COARSE 

MESH 


NUMERICAL 

REGION 


FINE 

MESH 


ANALYTICAL 


REGION 


*0(Wall) 


$  ■ 

(X,0) 


FIC.  1.  LATTICE  POINTS  IN  THE  PLANE 
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Subroutine  STEP  accepts  the  solutions  of  the  differ¬ 
ence  equations  from  the  main  program  at  horizontal  point  £n 
for  all  vertical  points  ij)^  and  then  steps  in  the  ^-direction 
to  calculate  the  species,  Y^,  and  the  enthalpy  ratio,  G,  at 
point  for  an  ^  points.  This  data  is  supplied  to  the  main 

program,  for  the  subsequent  calculation  of  the  other  thermo¬ 
dynamic  properties. 

The  mesh  in  tne  $  direction  (see  Fig.  1)  is  divided 
into  several  regions.  From  the  wall  at  if)  =  0,  c>n  analytical 
region  extends  to  $  =  if)l  in  which  species  and  enthalpy  ratios 
are  found  from  analytical  expressions.  For  the  region 
ipi  <  ij)  *  ip„,  the  difference  equations  are  solved  using  a 
fine  mesh  immediately  above  the  analytical  region  and  a 
course  mesh  in  the  remainder  of  the  region.  The  fine  mesh 
size  was  required  to  provide  the  necessary  numerical  accuracy 
in  the  solution  of  the  equations.  There  is  an  expansion  re¬ 
gion  above  $  =  $  . 

1.  Equations  for  Variables  Computed 
in  Analytical  Region 


Cl 


_  JLjsl 

3A  i 


1_ 


(6) 


m 


GAMMA  = 


(7) 


6 


BETA  =  (Cl) (P  )  {Gn+1  -  Gn  } 
e  l  o  o  J 


(8) 


ALPHA  = 


Gn+1] 
o  J 


-  a/0i  (GAMMA) 


(9) 


The  analytical  region,  extending  from  0  =  0 
(wall)  to  0  =  0l  may  be  subdivided  into  an  integral  number  of 
0  steps.  In  this  interval,  -.lie  enthalpy  ratios,  G^ ,  and 
species  mass  fractions  (Y,  ) ,  are  computed  from  the  relations: 


/. 


Gn'1’1  =  Gn+1  +  V01  \ALPHA  + 
l  o 


V 


BETA  (0.-0x )  _  ) 

- - - - -  +  y0i  GAMMA) 

-  (In  a)}  g.  j 


n+l_ 

i 


n+1 


1.0-{  i  (tn  a)}  g. 


(10) 


(11) 


The  static  enthalpy,  h,  is  found  from  the  stagna¬ 
tion  enthalpy  and  velocity  by: 


h.  =  h  G. 
l  e  l 


1 

2 


*2 

u .  u 

i 


2 

e 


(12) 
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The  mixture  temperature  is  then  obtained  from 
subroutine  ENTHLP  which  computes  temperature  from  static 
enthalpy  and  species  mass  fractions  (see  Ref.  2). 

The  mixture  molecular  weight  (WT)^  and  density 
ratio  (RH) ^  are  found  from 


(WT). 


1.0 


<Vi 


£ 

k 


“k 


(13) 


(RH). 


(WT)  i 


(14) 


<Ve 


th 

=  molecular  weight  of  k  specie 
=  species  reference  mass  fractions 


(at  edge) . 


and  (Y^)e  are  input  data  to  the  problem. 

The  g^'s  in  Eqs. (10)  and  (11)  and  the  velocity 
ratios  lu  are  computed  in  subroutine  HERB  using  the  equations 
of  Ref.  1. 

The  incompressible  viscosity  "I-VIS"  is  found 


The  quantities  “k 


from: 

(I-VIS) .  =  u./u.  . 

l  li 


(15) 
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If  the  compressible  viscosities  "C-VIS"  are  com¬ 
puted  (see  Sense  Switch  1  Option  in  Section  V) ,  they  are  found 
using  the  relation: 


(I-VIS) 
(C-VIS).  =  < - 


o 


€  Pp 

(M)i(p)i 


(16) 


The  parameter  (C/fi)  is  explained  in  Sections  II 


and  III. 

The  viscosity  p  and  the  parameter  |  are  computed 
from  the  relations: 


=  (3.05-lQ-e)T/a 
MAIR  T+111.0 


T  in  deg.  Kelvin 


(17) 


\2 


_ i 

1 . 0-C03  0 


si 

d 


<3X 


( In  a) 


(18) 


where  o  and  8  are  functions  computed  in  subroutine  HERB  (see 
Ref.  1) ,  and  the  subscript  s  which  is  explained  later  in  Section 
I  I A  refers  to  edge  of  the  sublayer. 


The  coefficients  of  these  equations  are  com¬ 
puted  and  then  the  resulting  set  of  linear  simultaneous  equa¬ 
tions  are  solved  in  subroutine  STEP,  first  for  species  and  then 
for  the  enthalpy  ratio.  The  incompressible  eddy  viscosities, 
u^,  are  supplied  to  STEP  by  subroutine  HERB.  These  viscosities 
are  then  corrected  at  each  $  point  for  compressibility: 


UCOMP  UINCOKP 


(20) 


where  u^  and  g^  are  also  obtained  from  HERB.  These  compressible 

u's  are  used  in  the  computation  of  the  a,  b,  c  coefficients  in 

the  generic  difference  Eq.  (19)  and  will  be  referred  to  as  u 

w- .h  no  subscript.  For  the  species  conservation  equation,  the 

P.  in  the  difference  equation  (19)  represents  (Y  )  and  the  coef- 
l  K  i 


ficients  are: 
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Se  (A  0) 


(21) 


ai  =  V%  +  Se(A  Wi  4  Un  0) 


(22) 


,  f .  ~n+l  ~n+ll 

bi  -  '  +  U1J} 


(23) 


'  Se(A*)<*)i.  dj  Un  CT) 


(24) 


<<V  i  -  *  V V  i 


(25) 


In  the  fine  mesh  region  (Fig.  1)  the  above 
relations  are  used  from  i  =  11  to  i  =  lr  with  A  ip  being  the  fine 
^-interval.  ipn  is  the  first  fine  mesh  point  above  ip1  and  $lr 
is  the  last  fine  mesh  point  below  ipa  . 

In  the  coarse  mesh  region,  the  above  rela¬ 
tions  are  used  from  i  =  3  to  i  =  M  with  (A  ip)  being  the  coarse 


tp  interval. 
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For  i{)a  , 


S  (A  0)s 
e  y  c 

2A  x 


(26) 


B  is  the  total  number  of  fine  mesh  intervals.  Equations  (22) -(25) 
become: 


aa 


~n+l 

U2-is 


+  f-  (A0)c  0a 


d?  Un  a) 


(27) 


ba  -  - 


xv  +  cfj  +  B 

Ya  2+H 


~n+l 
u  , 

2+h 


(28) 


c=  ■  -  f  (A  «c  *  4  (ta  a> 


(29) 


‘V.  =  -VV"  (30 

For  the  generic  form  of  the  enthalpy  ratio 

difference  equations,  relations  (21) -(24)  and  (26) -(29)  are 

used  "/ith  the  Schmidt  number  "S  "  replaced  by  the  Prandtl  num- 

e 

ber  ">  . " 
e 

Relation  (25)  for  the  right-hand  side  of  (19) 


is  replaced  by? 
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where 


a.  =  -  [*e<<  +  (R).  -  (R1..J 


(31) 


R.  = 

l 


u 


e  2h 


1  - 


e  i 


u .  , 

l+*5 


R 


-  u2 
i+2  l+l. 


(32) 


which  is  an  approximation  to  the  term  involving  — 


u 


liHl 


u 

e 

2h 

L  e 


1  - 


in  the  energy  equation  (1) .  For  the  Xg  in  Eq.  (31) 


Pe(A  *)2 
A  X 


(33) 


where  the  appropriate  A  $  is  used  depending  on  coarse  or  fine 
mesh  region.  For  i  =  2  Eq.  (31)  is  used  for  d2  with  X 

e 

replaced  by 


P  (A  0)2 

p  P 

2  (A  \) 


1.0 


+ 


1.0 

B 


(34) 


B  is  the  number  of  fine  mesh  intervals. 
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(b)  Boundary  Conditions  for  Difference 
Equations  at  i b  -  il) r 

At  0  =  the  generic  forms  of  the  species  and 
energy  equations  are  replaced  by  special  analytical  relations. 
For  the  species  equation,  these  are: 

ai s  =  0  (35) 


^is  " 


il) i  +  (A  0) 


u 


& 


2  (A  0  (A  ^)F  Se  A  I  2  (A  i) 


(36) 


$i  +  (  A  $) 


Un  0)]  [i-SCA-D  ] 


Cls  = 


~n+l 

u 


*4 


(37) 


+  ( A 

.  2  (A«)  . 


<vf + 


ii>i  +  (A  0) 


2 


t*-1-0]  ‘Vo 


+ 


1 

A 


3  [ 1~  ^  (tn  g)  'gQ 

SeV'2i?T  <p[^i  +  (A  $)p] 


N 

. 

/ 


(38) 
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where  (A  ip)„  is  the  fine  mesh  interval  and 


A  =  1.0  + 


[3  (A  £)]  [l.O  -  ~  (In  a)-gi] 


(39) 


After  solution  of  the  species  equations  for 
the  range  of  ip  from  0,  to  ip.,,  the  species  values  at  the  wall 
are  calculated  using  the  following  relation: 


(40) 


tine  ENTHLP  from 
temperature. 


The  wall  enthalpy  ratio  is  found  by  subrou- 

-  i  1 

the  wall  mass  fractions  (y  ),  and  the  wall 

O  K 


values, 


G^,  with 


The  energy  equation  is  then  solved  for  the 
the  following  boundary  conditions  at  ip  =  ipx  : 


aie  ~ 


0 


(41) 
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MA^)f  u?^1 
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4  72 
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3 „  Trea'.men  t  of  Psi  Expansion  Region 

To  ensure  that  the  solutions  satisfy  the  boundary 
conditions  at  the  edge,  an  expansion  region  is  included  in 
the  0  direction. 

Before  solution  of  the  species  equations  at  the 

current  step,  those  solutions  obtained  at  0=0  for  the 

LM 

previous  step,  are  compared  with  (Y.  )  ,  which  are  input.  If 
these  comparisons  differ  by  more  than  a  specified  tolerance, 
called  EPS,  for  any  one  of  the  species,  the  values  (Y.  )  .  are 

K  0 

prescribed  at  an  additional  ij)  point  which  is  added  to  the 
mesh  (TM  is  increased  by  one).  The  convergence  test  is: 


Is 


(Vur-<Ve 

<Ve 


*  EPS  ? 


(45) 


NO  -  add  1  point  to  mesh 

YES  -  do  not  add  a  point  to  mesh. 

A  similar  test  is  performed  on  the  energy  equa¬ 
tion  solution.  The  test  is  as  follows: 


I  s 


1.0 


*  EPS  ? 


(46) 


NO  -  add  1  point  to  mesh 


YES  -  do  not  add  a  point  to  mesh. 
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At  all  points  in  the  expansion  region,  the  u 

and  u  viscosity  parameters  are  set  equal  to  the  values  of  u 

and  u  at  )  =  0.  The  compressibility  correction  on  u  is  not 
M 

applied  in  this  region. 

4 .  Equations  for  Parameters  Computed  after 
Solution  of  Difference  Equations 

Upon  obtaining  the  solutions  to  the  species  and 

energy  equations  in  the  numerical  region,  the 

main  program  computes  mixture  temperature  ratios,  T./T  , 

1  a 

molecular  weights,  (WT) . ,  density  ratios,  (RH)^,  and  incom¬ 
pressible  viscosities,  (I-VIS)^,  using  Eqs.  (12)  through 
(15) .  If  desired,  compressible  viscosities  are  obtained, 
using  (16).  These  parameters  are  printed  as  output. 

The  program  then  computes  a  value  of  —  (in  !J)  to 
be  used  for  the  next  wtep  in  the  x  direction.  The  methods 
used  distinguish  the  substructure,  reference,  and  sublayer 
versions  and  will  be  detailed  in  Sections  II  and  III. 

Having  obtained  (^n  (7)  ,  the  physical  coordinate 
"X"  may  be  found  from  the  coordinate  X  using 


X-X. 

in. 


r* 
o  u  J 
°  e  *in. 


.3^  I  d 


i-ffl  |  (In  Q) 


dX 


(47) 


<7 

-  M, 
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where  X.  is  t  ^  initial  value  of  X,  and  x ■  is  the  corre- 
m.  Ain. 

sponding  value  of  X* 

Other  quantities  computed  and  printed  are  the 
heat  transfer  "Q-DOT,"  and  the  skin  friction  coefficient  CF, 
which  are  defined  by  the  following  relationships: 


.001285  \X 

0  1 
=  I 

p  u  h  a 

o 

FJ 

o  e  e 

J2  <P  P 

e 


(48) 


CF 


2.0  fi 

o  1 

o 

mJ 

<pa  T 


(49) 


The  units  of  q  are 


BTU 


ft  -sec 


Electron  mass  fractions  are  computed  from  the 


mass  fractions  of  Oa“  and  NO  using  the  following  relation; 


Y  _ 


(49a) 


The  electron  density  in  particles  per  cubic  centi¬ 
meter  is  computed  from: 

ELECTRON  PARTICLE  DENSITY  =  ( 5 . 67xlOa a ) ( p  _)(Y  _)  (49b) 

e  e 
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where  p^_  is  the  density  in  slugs/ft  . 

5 .  Finite  Rate  Chemistry  Option 

An  option  has  been  provided  in  the  program 
whereby  one-dimensional  finite  rate  chemistry  reactions  are 
computed  using  the  technique  of  G.  Moretti  (Refs.  3  and  4) . 
Using  this  option,  the  terms  containing  the  species  mass 
fractions  are  modified  to  reflect  the  coupling  of  the  one¬ 
dimensional  finite-rate  chemistry  relations,  with  the  two- 
dimensional  diffusion  equations.  However,  since  the  correct 
time  step  for  the  chemistry  equations  is  not  known  until  the 
diffusion  equations  are  solved,  a  time  step  iteration  is 
required. 

The  iteration  procedure  is  as  follows: 

(1)  An  approximate  A  X  is  c  imputed  from  the 
previous  temperature  and  species  profiles: 


(A  X) 


(x)  _ 


P  u 
o  e 


i-<p3e 


LdX 


(£n  a) 


a 


(A  X) 


(50) 


(A  t) 


0 )  _  (A  X) 


(i) 


(P  u  //*  )u. 

e  e  e  i 


(51) 


20 


(2)  A  finite  rate  chemistry  step  is  performed  at 
each  mesh  ooint,  $  ,  on  the  species  mass  fractions,  using  the 
corresponding  temperature  and  density  profiles  and  the  "At" 
profile  computed  from  Eq.  (51). 

(3)  The  diffusion  equations  are  then  solved  using 
the  chemically  modified  species  profiles. 


(4)  New  values  of  *r~  (tn  a)  and  c/u  are  computed, 

dX 
(a) 


and  then  (AX)  found,  using  (50). 


(5)  (A  X)  ^  is  compared  to  (AX) 


If  they 


are  within  the  specified  tolerance,  the  species  and  energy 
solutions  are  printed  as  the  correct  solutions.  If  the  toler¬ 
ance  is  not  met,  steps  2-5  are  repeated  with  the  new  value  of 
(A  t)  .  Note  that  in  this  case  the  new  value  of  —  (£n  cj)  from 
step  (4)  is  not  used  in  the  left  side  matrices  of  the  diffusion 
equations  for  the  next  iteration.  It  is  only  used  to  get  a  new 
(A  t)  approximation  for  the  chemistry  calculation. 

B.  Numerical  Methods  of  Solution  of 
Basic  Equations 

Partial  differential  Eqs.  (1)  and  (2)  are  solved  by 
an  implicit  second-order  central  difference  method  known  as 


the  Crank-Nicolson  Difference  Equation. 
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FIG.  2.  CRANK-NICOLSON  LATTICE  POINTS 


Assuming  a  two-dimensional  mesh  of  lattice  points 
with  "n"  representing  the  horizontal  or  x  direction  and  "i" 
representing  the  vertical,  or  0  direction,  we  solve  the 
equation, 

dp  d  I*-  dpi  . 

dx  “  W  |_U*X'^  dtbj  ( 

at  a  point  (n+1) ,  i,  assuming  known  values  of  P  at  point  n 
for  all  values  of  i.  Replace  the  right  side  of  (52)  by  a 
linear  relation  for  the  derivative  from  i-^  to  i+^: 


Each  term  in  the  numerator  of  (53)  is  approximated  in  the 


same  manner. 


1  -s  _ 

n+1  n+1 

~P.  -P.' 

[  ~  dp 

l+l  l 

r 

=  u 

i+%  i+i2 

A  0. 
i 

n+1 


u 


dP 
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n+1 


i-k 


~n+l. 
-  u 


p.-p.  . 

i  i-l 


A  0. 


n+1 


The  left  side  of  (2)  is 


dP 

ax 


n+1  n 
P.  ~P. 

l  l 


A  X 


Inserting  (54) -(56)  into  (53) 


n+1  n 
P.  -P. 

l  l 


A  X 


(54) 


(55) 


(56) 


(A  0)T  !  ui+^  (pi+rpi} 


i  L 


u.  ,  (p.~p.  J 

x-H  x  i-l 


n+1 


•  (57) 


Multiplying  both  sides  of  (57)  by  (A  02 )  and  rearranging  terms 


the  difference  equation  becomes: 
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I 

i 

I 

I 

i 

I 

I 

I 

I 

I 

I 

« 


~  n+l 

U.  L,  P.  .■ 

i-h  i-i 


(A  4>)3, 

JL  *** 

— 7 -  +  u  •  .i  +  u  ■  L 

A  \  i+h  i -H 


n+r  ~  n+l 
P.  +  u.  .  P . . . 

l  i+*5  l+l 


<*»>!  n 

i  ‘ 


Let 


(58) 


a.  =  u 


i-*a 


b. 

l 


(A  *)a 

2_ 

- 1 -  +  u.  .  +  u.  , 

A  X  i+J5  i-*5 


c .  =  u . 


i  i+h 


d. 

l 


^A  ty)2. 

1  pn . 


Ax  i 


Then  (58)  may  be  written  in  the  form. 


n+l  n+l  ,  _n+l  , 

a.P.  .  +b.P.  +  c , P ,  .  =d.  . 

i  i-i  ii  i  l+l  i 


(59) 


Since  all  for  a  particular  mesh  line,  n+l,  are  solved 
simultaneously  for  i=l  to  l,  then  we  have  a  set  of  l  equa¬ 
tions  of  type  (59)  for  the  unknown  P^‘s. 
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written  in  matrix  form  as 


bx  cx  0  0  0  0  0 

i 

•-» 

aa  bg  ca  0  0  0  0 

p3 

da 

0  a3  b3  c3  0  0  0 

•  •  t  •  *  •  ■ 

0  0  0  0  0  a.  b. 

P, 

d3 

L.  l  i_ 

_  1  - 

t  _ 

(61) 

or 

AP  =  D  (62) 


To  solve  for  the  unknown  P's,  the  coefficient  matrix 
(called  A)  is  factored  into  a  product  of  two  matrices  as 
follows: 


[M  N]  P  -  D 


(63) 


-mm 

I 

f 
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where 


M  = 


&  0  0  0  0 

Os  &  0  0  0 

0  Gta  0  C 


(64) 


1  Vi  0  0  0 

oi  ya  o  o 

oo  l  y3  o 


(65) 


L  1  J 

The  0L,  /3,  and  >  values  of  M  and  N  can  be  evaluated 
by  multiplying  M  and  N  and  setting  the  elements  of  this  product 
matrix  equal  to  the  corresponding  elements  of  A.  V?hen  this  is 


done  it  is  found  that 
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and  =  ax  =  0,  /3X  =  b1 ,  y1  =  c  ,  /t'i  . 

In  (63)  let  Y  =  NP.  Then  since  M  is  a  bi-diagonal 
matrix  the  transformed  unknown  column  matrix  Y  can  be  solved 
recursively  from  j=l  to  l,  as  follows: 


IKWt 
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yi  =  dj  /  $i 

ya  =  (da-ttay:)/^  (67) 

yt  =  • 


At  this  point,  a  boundary  condition  is  imposed  and  y^  is 
modified  such  that 


,  y* 
l+*; 


(68) 


where  c,  =  u  . 
t  e 

The  solutions  P.  may  now  be  calculated  from  y  by 

1.  **■ 

sweeping  backward  from  t  to  1 


28 


y 


i 


p, 

i 


yi-Vi+i' 


i=t-l, 1-2, 


1. 


I 

1 

I 

I 


I I .  SUBSTRUCTURE  AND  REFERENCE  HYPOTHESES 

A.  Calculation  of  d/dx(^n  ff) 

Integrals  are  computed  over  the  variable  "ZETA" 

(which  is  a  transformed  ifr  coordinate)  fo  temperature  T, 

species  Y  ,  and  normal  coordinate  "YCORD,,"  These  integrals 
K 

are: 


-  jl  r430  . 

Ts  430  J  Tzd^ 
o 


(70) 


‘Vs  =  430  J  <V{  dC 

o  b 


(71) 


where  subscript  s  denotes  some  mean  value  across  the  layer, 
As  proposed  by  Coles,  the  substructure  hypothesis 


(72) 


where  /ig  is  the  mean  value  of  viscosity  in  the  region 
0  s  C  ^  430  and  /i  is  the  incompressible  viscosity  inde¬ 
pendent  of  the  variable 


-4 
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The  normal  y-coordinate  is  obtained  using: 


YCORD  = 


1  a  J 

=  P  M 
H  ere 


e 

-  dz 

? 


(73) 


The  computing  scheme  is  then  as  follows:  The  finite 


difference  equations  have  been  solved  for  the  properties  at 
Xn+1  using  the  value  of  ~  (£n  (7)  at  X*1  (see  Section  I,  B)  . 

'•*  A 

(For  the  first  step,  —  (£n  a)  is  an  input  value  to  the  pro- 

Q  A 

gram.)  Having  now  the  vahae  of  from  Eq.  (17),  the  value 

~  (tn  cr)  to  be  used  for  the  step  XA+1  to  xn+2  is: 


(In  cr)  =  -  — 


n+1  n 
s  s 


n  I  n-f  l  n 

M  Lx  -x 


(74) 


The  value  of  /x  for  n  =  0  is  computed  from  Eq.  (17)  using  the 

s 

input  temperature  and  species  profiles. 


(tn  (7)  thus  lags  the  remainder  of  the  solution  by 


one  step. 


Referring  to  Fig.  1  of  Section  I,  the  fine  mesh  region 
for  the  substructure  hypothesis  extends  over  the  interval  ($x,$2) 
The  p  step  size  between  Pi  and  ips  is  therefore 
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4> 

(A  ip)  =  — — —  (75) 

F  tr 

A  \ 

where  K  is  an  input  number  to  the  program. 

The  integrals  over  the  logarithmic  region  [Eqs.  (70) 
and  (71)]  are  found  by  trapezoidal  quadrature  in  subroutine 
INTEG,  over  the  values  of  £  used  in  the  finite  difference  mesh. 
Temperature  and  species  values  for  the  upper  limit  C  =  430  are 
found  by  linear  interpolation. 

Since  the  definition  of  £  changes  at  £  -  10.6  (see 
Appendix  III,  Ref.  1),  a  special  approximation  scheme  was 
used  for  the  £  interval  bracketing  10.6.  This  interval  was 
split  into  two  intervals  namely,  Cl  to  10.6,  and  10.6  to  £a  , 
where  Ci  and  Ca  are  the  mesh  values  of  £  that  bracket  C  =  10.6 
(see  Fig.  3)  . 


FIG.  3.  QUADRATURE  DIAGRAM  AT  C  =  10-6 
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A  step  drop  in  temperature  is  imposed  at  £  =  10.6,  yielding 

two  temperature  values,  denoted  T*  and  T '  _ .  Two  trape- 

10 . b  10 . b 

zoidal  integrations  are  then  performed,  the  first  from  Ti  to 
+  _ 

T1r>  c/  the  second  from  T1(.  to  Ta  .  (See  shaded  areas  I  and 
10 . b  10 . b 

II  in  Fig.  3 . ) 

The  values  of  T+  -  and  T,_  _  are  obtained  as  follows 

10.6  10.6 

T10.6  =a+b^+c^S  (76) 

where 


0  =  ip/10.6  (77) 

A  =  T  (78) 

o 


C  = 


rZi. . 

-Ui 

=S_  +  T  { 

ua  o  l 

_1_  _ 
ua 

ui  -  u2 

Ta-To 

B  : 

u3 

C  ua 

T“ 

10.6 

is  found 

from 

a  backward 

(79) 


(80) 


extrapolation  of  the  temperatures  Ts  and  T^  ,  where  T43Q  is 


the  temperature  previously  found  by  interpolation  at 


c  =  430. 
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B .  Modification  of  Grid  Mesh 

in  Normal  Direction 

The  number  of  grid  mesh  points  in  the  normal  or  0 
direction  is  determined  by  the  value  of  0  ,  or  upper  limit 
of  0  in  the  numerical  region  (see  Section  I) .  The  initial 
value  of  0  is  known  and  a  0  spacing  of  di  /N  is  input  as 
(A  ^)c*  As  calculation  proceeds  in  \  direction,  in¬ 
creases  and  additional  mesh  points  are  added  with  the  (A  $)c* 
When  the  total  number  of  these  points  reaches  2N,  the  pro¬ 
gram  automatically  doubles  (A  t)  and  halves  the  number  of 
points,  keeping  solution  values  at  every  alternate  point  of 
the  original  A  0. 

For  the  fine  mesh  region,  whose  interval  is  also 
doubled,  alternate  points  are  kept  for  the  lower  half  of  the 
region.  Points  for  the  upper  half  of  the  new  fine  mesh  region 
are  linearly  interpolated  (see  Fig.  4) . 

For  the  coarse  mesh  region,  every  other  point  of  the 
original  mesh  is  retained. 
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C.  Treatment  of  Wall  Chemical  Reactions 

Finite  rate  chemistry  reactions  are  computed  at  all 
iji  points  except  at  the  wall,  where  an  equilibrium  chemistry 
computation  is  performed. 


Reference  Method  Option 


When  the  reference  method  option  is  exercised, 
cl 

(^ii  or)  is  set  equal  to  zero.  The  reference  state  is  still 
given  by  the  mean  substructure  values,  i.e.,  T  ,  (Y,  )  ,  , 

Cf/jl,  and  YCORD  using  relations  (70),  (71),  (17),  (72),  and 


(73)  respectively. 
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III.  SUBLAYER  HYPOTHESIS 

A.  Calculation  of  d/dXC^n  g) 

The  sublayer  assumption,  proposed  by  Baronti  and 
Libby,  asserts  that  the  Reynolds  number  based  on  the  height 
of  the  laminar  sublayer  is  an  invariant  of  the  compressibility 
transformation.  Instead  of  relation  (72)  for  cr/^t,  one  uses 
the  following: 


a  1  r.56*1.8  Pg 

f  =  loTe  JQ  T  72? 


(81) 


and 


H  .56.18 
s 


10 


dX 


(£n  <7)  = 


s_  r  d_  rs 

.6  J  dX  P 

o  ' 


p-'  .  a.  ,  _l_  r56-18 


M 


dX 


(Ms}  10.6  I 
o 


M 


(82) 


where 


_d_ 

dX 


/  p  1 

P 

fP 

_S  1 

_ S 

d 

e 

Ip  1 

P 

dX 

P 

e 

_s  _d_ 
P  dX 


IP  1 

e  1 

P  j 

s' 

— 

(83) 


d 

|Pe| 

1  ~  — - 

'Pe| 

dT 

+  V 

3 

r  P 

e 

dX 

P 

1  3T 

P 

dX 

L 

k 

SYk 

i  p 

d_ 

fPe’ 

_  _B_ 

lPe 

dT 

+  y 

d 

f  Pe  | 

dX 

lpsi 

St 

>Ps 

dX 

L 

k 

3Yk 

dY 

} 

dX 


dY 

J 

dX 


(84) 


(85) 
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dT 


I  <■ W  s 

JL_  k _ 

’s  l  We 

k 


(96) 


_L  lie]  =  _J: _ 1 _ 

VWTe  (Mk>s  1  (Ws 


(87) 


d_ 

dX 


Idhl 

f  dT 

dT  j 

idX  , 

(88) 


3.05x10  8t1 

L  ,  .  T 

T  +  111.0  j 

\J"~’  ~  T  +  111.0 

(89) 


The  derivatives  dT/dx  and  dY^/dX  are  evaluated  numerically 
as : 


dT  =  (T) n+1- (T^ ” 

dx  •  xn+i-xn 

1  —  X  1 

dX 


<V  -(Yk> 

xn+1-xn 


(90) 


(91) 


In  relations  (81)  to  (91) ,  the  subscript  s  refers  to 


parameters  evaluated  at  0  =  56.18. 
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The  sublayer  hypothesis  applies  to  the  fine  and 
coarse  mesh  ip  regions,  just  as  the  substructure  hypothesis 
does.  However,  the  fine  mesh  region  now  extends  from  ipi  to 
Ip  =  56.18  and  can  be  divided  into  a  prescribed  number  of 
•>  ntervals . 

The  normal  y-coordinate  is  still  computed  from 
relation  (73).  As  in  the  substructure  hypothesis,  calcula¬ 
tion  of  d/dx(^n  (7)  lags  the  remainder  of  the  solution  by 
one  step. 

B.  Modification  of  Grid  Mesh 

in  Normal  Direction 

When  the  criterion  for  halving  the  number  of 
intervals  in  the  ip  direction  is  inet  (see  paragraph  1, 

Section  II,  B)  ,  the  fine  mesh  ip  points  of  the  original  solu¬ 
tions  are  retained,  i.e.  the  fine  mesh  region  is  undisturbed 

for  0  *  0  £  56.18  (ip  =  56.18). 

s 

For  the  coarse  ip  mesh  region,  56.18  <  ip  £  ip  ,  every 
other  mesh  point  of  the  original  solution  is  retained. 

C.  Treatment  of  Wall  Chemical  Reactions 

Finite  rate  chemistry  reactions  are  computed  at  all 
Ip  points.  At  the  wall,  a  t^e  step  is  used  that  is  equal  to 
uhe  time  step  computed  at  the  ip  value  closest  to  the  wall. 
See  Eq.  (51)  for  definition  of  time  step  (A  t)^. 
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IV.  DESCRIPTION  OF  INPUTS 


A .  Calculation  of  Initial  Input  Data 
1 .  Given  Information 

The  following  information  must  be  specified: 

a.  External  conditions  (u  ,  p  ,  u  ,  T  ,  Y,  ) 

e  e  e  e  ke 

k  =  .1,  ...  7  representing  species  0,N,N0,  Oa  ,  Oa  .  N2  ,  and  N0+. 

b.  Wall  conditions  (T  ,  (Y.  )  )  k  =  1,  ...  7. 

o  k  o 

(1)  Initial  compressible  skin  friction 

2 

coefficient  c..  =  T  /p  u  . 

f  o  e  e 

(2)  Initial  compressible  Reynolds  number 
based  on  momentum  thickness  - 


Rfl  H  P  u  • 

0  e  e  e 

c.  Initial  temperature  variation  with  velocity 
ratio  T(u/ue)  through  viscous  layer. 

d.  Initial  species  mass  fraction  variation  with 

velocity  ratio  Y,  (u/u  )  k  -  1,  ...  7  through  viscous  layer. 

K  0 


2 .  calculation  of  (7//i 

As  a  first  step  in  determining  the  input  data  the 
parameter  cr/^i  must  be  related  to  the  incompressible  skin  fric¬ 
tion  coefficient  Cf  .  This  procedure  varies  depending  on 


whether  the  substructure  or  sublayer  hypothesis  is  utilized. 
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a .  Calculation  of  CT/ji  According 
to  Substructure  Hypothesis 

For  the  substructure  hypothesis,  take 


*  =  4  =  h 
c  *s  p 


T  (Y,  ) 
s  k  s 


(92) 


where  ^  denotes  the  viscosity  of  the  mixture  evaluated  at 
s 

the  temperature  T  ,  and  composition  (Y,  )  where  these  latter 

S  K  S 

are  given  by 

,  .430 

Ts  =  j-  j  TdC  (93) 

o 


<Vs 


i  r 

430  J 
o 


430 


YkdC 


(94) 


k=: 


.7 


In  accordance  with  items  c  and  d,  T  and  Yn  are  known  func- 

k 

tions  of  u/ue>  Furthermore,  from  the  Eqs.  (AIII-1)  through 
(AIII-6)  -given  in  Ref.  1  and  the  relations 


<p  = 


(95) 


and 


C6  “  exp  +  2'03 


(96) 


u  1  u 

u  u 

e  T 


(97) 
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one  can  obtain  the  correspondence  between  the  velocity  ratio 
u/u  and  C  for  any  particular  value  of  C,.; 

6  X 


c 


(98) 


Tnus  the  integrals  appearing  in  (93)  and  (94)  can  be  evaluated 
(numerically  if  necessary)  and  will  depend  only  on  a  choice  of 
Cf*  Thus  also  Cr /\i  will  be  related  uniquely  to  C^;  i.e.  : 

ff  =  Cf  (Cf)  .  (99) 


b.  o/li  According  to  Sublayer 
Hypothesis 

For  the  sublayer  hypothesis,  take 

r6-18ls 

77  10.6  J  P 

*“  o 


(100) 


where  p  and  fi  denote  the  density  and  viscosity  of  the 
s  s 

mixture  at  ip  =  56.18.  In  general,  the  density  is  related  to 
the  temperature  and  species  by 


P 


(101) 
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where  the  denote  the  molecular  weights  of  the  individual 
species  and  are  given.  Since  the  and  T  are  related  to  the 
velocity  ratio  through  c  and  d  above  we  have 


P 

P 


Nov/  relate  PQ/P  to  0  by  the  relation 


(102) 


(103) 


so  that,  as  in  the  sub-structure  case,  there  can  be  written 
formally 

a  cr  — 

=  =  Z  (Cj  •  (104) 

3 .  Calculation  of  <p 

a.  With  Cf  given,  solve  for  Cf  (by  iteration) 
from  the  following  equation 

ff  =W£  _ 

Cf  Pe  M 


(105) 
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b.  With  Ra  given,  use  Eq.  (105)  and 


R  n  i  — * 

8  _  1  g 

R  "  M 
7T  e 


(106) 


Cf  =  Cf(R  ) 


(107) 


where  Eq.  (107)  is  given  graphically  in  Fig.  5.  The  solu¬ 
tion  for  Cf  is  again  obtained  by  iteration. 

4.  Calculation  of  <p,  ij) 

_ 0 _ M 

Once  an  initial  value  of  has  been  obtained 
the  corresponding  values  of  <p  and  Cg  follow  from  (95)  and 

(96),  while  is  obtained  from 

M 

=  -  30.81  +  2.43  Ci  Ci  +  2.47  Cl  +  (C§”Ci  )<P  +  1-7  Cg 

(108) 


where  Ci  =  0.131  Cg« 


5 .  Calculation  of  Initial  Profiles 
T  W)  ,  Yk(ft),  G(VQ 

From  the  given  inputs  there  is  available  T(u/u^) 
and  Y,  (u/u  ).  G(u/u  )  is  obtained  from 

K  G  6 
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G  = 


H 


u 

u 


2  U 


e 

2H 


From  the  parametric  relations 


—  =  =-  =  function  of  £ 
u  u 
e  e 


0  =  0(0 


I 


I 


I 

1 


given  in  Appendix  III  of  Ref.  1  there  is  tabulated  the  rela¬ 
tion 

u/u  =  u/u  (0) 
e  e 

where  G,  T  and  Y,  can  be  obtained  as  functions  of  0.  These 

K 

are  plotted  in  graphical  form  from  which  the  desired  values 
corresponding  to  the  previously  selected  0-mesh  points  0^ 
ere  read  off. 

6 .  Numerical  Example  (Sublayer  Hypothesis) 
a.  External  conditions 

p  =  1.344  x  10  3  sluqs/ft3 
e 

u  =  2120  ft/sec 
e 

T  =  12  2°  K 
e 

u  =  .4926  x  1.153  x  1C  8#/ft-sec 
e 


nr 
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Combining  e 


(YB)e  =  .232 
(Y-  )  =  0.768 


y  undissociated 
air 


(Y.  ) 

K  e 


=  0;  k  =  1,2, 3, 4,7  / 


b.  Wall  conditions 


T  =  306  K 
o 


(Y,  )  =  (Y,  )  ;  k  -  1,  ...  7  . 

ik  O  J\  6 


c.  Skin  friction  coefficient 


Cf  =  .0013  . 


d.  Temperature  distribution 


T=T  +  (T  -T  )  —  -  (T  -T  ) 

o  s  o  u  s  e  u 

e  e  e  c 


(Crocco  integral;  Pg  =  1) 


e.  Species  distribution 


Y,  (u/u  )  =  constant  =  (Y.  ) 

K  6  K  6 


,  (101)  and  (103)  gives,  using  the  numerical  data 
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2.51  +  0.22 


2.51  + 


2.33 

O 


194 

<ps 


so  that  from  (100) 


10.6 


2.51 


2.51 


1.165  64.67 

+  0  ~  <P2 

,  2.33  _  194 
0  ~  02 


Now  take  an  initial  guess  of  Cf 


2x10  3  so  that  from  (95) 


for  which 


0  =  22.4 


P 

~  =  0.451 


2.06 
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ind 


~=~-  =  1.095 
M 


From  the  given  data 


P  M 
o  o 

P  M 
e  e 


=  1.131  so  that  from  (105) 


PM  M  M  a 
o  o  e  's  — 

Cf  "  P  M  M  M  Cf 
e  e  s 


.0012  <  .0013  . 


A  second  guess  of  Cf  =  .0025  yields 

Cf  =  .00163  >  .0013  . 

A  linear  interpolation  gives  as  a  third  guess  Cf  =  .00212 
for  which  one  obtains 

Cf  =  .00131  «  .0013 

which  is  the  required  result. 

B.  Input  Formats  for  IBM  Programs 

In  this  section,  the  input  formats  for  each  of  the 
two  program  decks  -  Substructure  /  Reference  Hypothesis,  and 
Sublayer  Hypothesis,  will  be  described  in  detail.  Refer  to 
the  section  on  Nomenclature  for  allied  information. 
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The  terra  "card"  refers  to  the  standard  IBM  data 
processing  card  consisting  of  12  rows  and  80  columns.  The 
terra  "format"  refers  to  the  mode  of  input.  Symbolically, 
these  modes  may  be  defined  as  follows: 

I  integer  ±  XX  (no  decimal  point) 

E  floating  ±  X.XXX  ±  YY  (YY  is  the  exponent 
point  to  the  base 

10.  ±  X.XXX. 10 (YY) ) . 

For  the  E  mode,  the  decimal  point  may  be  shifted  from 
the  position  indicated  in  the  above  example  and  the  maximum 
number  of  significant  figures  is  governed  by  the  field  width 
assigned  for  each  "word"  of  data.  The  plus  (+)  sign  may  be 
omitted  in  all  cases,  except  for  the  sign  immediately  pre¬ 
ceding  the  exponent  for  the  E  mode.  An  additional  format  is 
the  Hollerith  mode  which  consists  of  alpha-numerical  informa¬ 
tion,  and  for  our  purposes,  will  be  utilized  exclusively  for 
an  identif icr t ion  input  card,  which  wilj  subsequently  be 
printed  as  a  title  at  he  head  of  the  output  listing.  It  is 
good  practice  to  "right-adjust"  data  words  within  the  indicated 
field;  that  is,  the  word  must  be  ehifted  to  the  extreme  right 


of  the  field. 
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1 .  Substructure  and  Reference  Hypotheses 


CARD  COLUMN 


DESCRIPTION 


1  1  Punch  the  number  0 


FORMAT 


2-72  Title  information  H 

2  6  Punch  the  number  "1" 

2  14-15  M  -  Number  of  coarse  mesh  points  in  "0"  I 

direction,  *  40 

Total  number  of  species  =7  I 

55  Punch  the  number  "1"  I 

60  Maximum  number  of  iterations  on  "Ax" if  I 

finite  rate  chemistry  option  is 
requested  5) [See  Eq.  (50)] 

64-65  Number  of  fine  mesh  points  in  "0"  I 

direction,  ^  25 

69-70  m  -  Print  cycle  nrunber  -  print  properties  I 

at  every  x_step,  5  10 

3  4-5  M  -  Punch  same  number  as  in  cols.  14-15,  I 

card  2 

10  Chemistry  option:  Punch  "1"  if  finite-  I 

rate  chemistry  requested;  punch  "0"  if 
no  chemistry 

15  "G"  input  profile  option  I 

"1"  -  input  temper-  inputs  on 


ature  in  °K  > 8th  set  of  cards 

"2"  -  input  stagnation  I  A+l  through  C 
enthalpy  ratios  / 

20  Reference  or  substructure  hypothesis  option  I 

"0"  -  substructure  hypothesis 
"1"  -  reference  hypothesis 
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CARD 

COLUMN 

DESCRIPTION 

FORMAT 

4 

1-15 

"A  0" 

for  coarse  mesh  =  (0M-0X )/M 

E 

31-45 

Final 

value  of  "X" 

E 

61-75 

"  _ 

o 

initial  <fi 

E 

5 

1-15 

"(C6)o 

"  -  initial  £6 

E 

.16-30 

"HE"  - 

reference  stagnation  enthalpy 
(ft3/sec3 ) 

E 

31-45 

"A 

-  CSI  step  size 

E 

46-60 

"DPSY" 

-  "0"  step  size  between  wall 
and  "0X " 

E 

6 

1-15 

"TROLL 

"  -  tolerance  for  iteration  on 
"A  X"  for  chemistry  (approx. 
0.05)  [see  Eq.  (50)] 

E 

16-30 

"EPS" 

-  tolerance  for  adding  point  to 
species  or  energy  solution 

matrix  (approx.  0.001)  [see  Eqs. 
(45)  and  (46)] 

E 

7 

1-15 

"  p  "  _ 

e 

Prandtl  number 

E 

16-30 

"S  "  - 
e 

Schmidt  number 

E 

31-45 

"V  " 

first  numerical  0  value  after  wall 

£ 

46-60 

nip  » 

e 

reference  temperature  (°K) 

E 

61-75 

"u  "  - 
e 

reference  velocity  (ft/sec) 

E 

8 

1-15 

n  p  "  _ 

e 

reference  density  (slugs/ft3) 

E 

16-30 

"  M  "  _ 

e 

reference  viscosity  (lb-sec/ft  ) 

E 

46-60 

d 

Initial  value  of  ^  (£n  (J)  (sub¬ 
structure  version  only) .  Use  zero  it  not 
krown. 

E 

73-75 

"SS"  - 
(col  s. 

punch  1.0  if  ^  (In  a) 

50-60)  is  not  zero 

E 
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CARD  COLUMN 


DESCRIPTION 


FORMAT 


9,  .  .A 


1-15,16-30, 

..61-75 
1-15, 16-30 


Initial  wall  species  for  species 
1  to  7 .  The  species  are  0,N, 

NO, Ob, Os, Ng  and  NO  • 


E 


A+1,A+2..B  1-15,16-30, 

..61-75 

1-15,  etc. 


Values  of  Yl  at  all  fine  mean 
points  along  initial  $- mesh 
line  from  $=$ x  to  point 
immediately  below 


E 


B+l, B+2 . . C  1-15,16-30, 
..61-75 

1-15,  etc. 


Values  of  Yj  at  all  coarse  mesh 
points  along  initial  ^-rnesh 
line  from  $=0a  to 
inclusive 


E 


Repeat  cards  "A+l"  to  "C"  for  remaining  species 
Ya  through  Y7 ,  and  then  for  "G"  profile.  Thus 
there  are  8  sets  of  cards  designated  A+l 
through  C 


C+l. .D 


D+l 


D+2. .E 


1-15,16-30, 

..61-75 

1-15,16-30 


(Yk)e  for  k  =  1, 2, . . . , 7 
(Mass  fractions  at  edge  of 
boundary  layer) 


E 


Wall  temperature  function  vs.  "X"  E 


(Tw 

)l  =  + 

Bx 

X, 

for 

X 

(Tw 

)8  -  Aa  + 

Ba 

X, 

for 

X 

*  Xi 

1- 

-15 

Ai  , 

(°K) 

16- 

-30 

Ag  , 

(°K) 

31- 

•45 

Bi  , 

°K/units 

of 

Xi 

46- 

-66 

Ba  , 

°K/units 

of 

Xi 

61- 

-75 

xiv 

Peue/ Me) ' 

where 

Xi 

is 

in  feet 

l-rl5, 16-30,  Molecular  weights,  M^,  for  species  E 
..61-75  1  to. 7 


1-15,16-30  These  values  are:  Mi =16.0,  M3sa14.0, 
M»— 30.0,  M4“32.0,  Mr  =32.0,  Me-28.0, 
M?  =30.0 
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CARD 

1 

2 


2 .  Sublayer  Hypothesis 


COLUMN 

DESCRIPTION 

FORMAT 

1 

Zero 

2-72 

Title  information 

H 

4-5 

Number  of  fine  mesh  points  in 
direction,  s  25 

I 

9-10 

M  -  Number  of  coarse  mesh  points  in  "0" 
direction,  s  40 

I 

14-15 

TOTAL  number  of  mesh  points  in  "0" 
direction  (sum  of  values  in  cols.  4-5 
and  9-10) 

I 

19-20 

M  (same  as  in  cols.  9-10) 

I 

24-25 

Punch  same  as  in  cols.  14-15 

I 

30 

"G"  input  profile  option 

Punch  "1"  -  if  temperature  in  degrees 

Kelvin  are  input  on  8th  set  of 
cards  A+l  through  c 

Punch  "2"  -  if  stagnation  enthalpy  ratio 

(h/h  )  are  input  on  8th  set  of 
cards  A+l  through  C 

I 

35 

Punch  a  "1" 

I 

40 

Punch  a  M7" 

I 

44-45 

m  -  Print  cycle  number  -  print  properties 
at  every  mth  x_step,  s  10 

I 

50 

Chemistry  option:  Punch  "1"  if  finite-rate  I 

chemistry 

requested; 

Punch  "0"  if  no  chemistry 
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CARD 

COLUMN 

DESCRIPTION 

FORMAT 

2 

55 

Punch  "1" 

60 

Chemistry  option  for  (OgT) 

"4"  -  (Oa)  is  not  included  in  chemistry 
reactions 

"5"  -  (0a)  is  included  in  chemistry 
reactions 

I 

3 

11 

Punch  "1" 

I 

20 

Maximum  number  of  iterations  on 
"A  X"  if  chemistry  option  is 
requested,  *  5  [see  Eq.  (50)] 

I 

4 

1-15 

"A  ip"  for  coarse  mesh  region  between 

Ip  =  56.18  and  Ip  =  y>M 

E 

31-45 

"£p"  =  final  value  of  CSI 

E 

61-75 

"<p  "  -  initial  <P 

E 

5 

1-15 

"C60"  “  initial  ZETA  DELTA 

E 

16-30 

h  -  reference  stagnation  enthalpy 
(ft*/see3) 

E 

31-45 

"A  £"  -  £  step  size 

E 

46-60 

"DPSY"  -  A  ip  ior  "ip"  between  wall  and 
ip  =  ipi 

E 

6 

1-15 

"TROLL"  -  tolerance  for  iteration  on 

"A  X"  if  chemistry  option  is 

requested,  approx.  .05  [see 

Eq.  (50)] 

E 

7 

16-30 

1-15 

"ERS"  -  tolerance  on  D  for  adding  pt.  to 
solution  matrix,  approx.  0.01  [see  Eqs. 
(45)  and  (46)] 

P„  -  Prandtl  number 
e 

E 

E 

16-30 

Se  -  Schmidt  number 

E 

31-45 

ipl  -  First  numerical  " ip "  value  above  wall 

E 
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CARD 

COLUMN 

DESCRIPTION 

FORMAT 

7 

46-60 

Tg  -  reference  temperature  (°K) 

E 

51-75 

u  -  reference  velocity  (ft/sec) 

E 

8 

1-15 

p  -  reference  density 
(slugs/ft3) 

E 

16-  30 

(X  -  reference  viscosity 
(lb-sec/fts) 

E 

46-60 

Initial  value  of  'In  0) 

QX 

?use  zero  if  not  known) 

T‘ 

61-75 

"SS"  -  punch  "1.0"  if  cols.  46-60 
is  not  zero  or  blank 

E 

9,  .  .A 

1-15,16-30., 

. .61-75 

1-15, 16-30 

Initial  wall  species  for  species 

1  to  7 .  The  species  are  ON  NO. 
Oa"  Oa  •  N3  and  N0+ 

E 

A+1.-A+2  ,  .B 

1-15,16-30 
. .61-75 

1-15,  etc. 

Values  of  Yi  at  all  fine  mesh 
points  along  initial  0-mesh  line 
from  0=0i  to  point  immediately 
below  0=56.18 

E 

B+lr  B+2 . .C 

1-15,16-30, 

.  61-75 

1-15,  etc. 

Values  of  Yi  at  all  coarse  mesh 
points  along  initial  0-mesh 
line  from  0=56.18  through  0=0M 

E 

Repeat  cards  A+l  through  C  for  remainirg  species 
Ya  through  Y?  and  then  for  G  profile  thus  there 
are  8  sets  of  cards  designated  A+l  through  c 


(Yk)e  for  k  =  1. 2  .  . 7 
(Mass  fractions  for  edge 
of  boundary  layer) 


C+l  D  1-15,16-30 

. . 61-75 

1-15  16-30 
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CARD 

COLUMN 

DESCRIPTION 

FORMAT 

D+l 

Wall  temperature  function  vs.  "X" 
(Twh  "A*  +  Bj.  X,  for  X  *  Xi 

(TV) a  =  Aa  +  Ba  X,  for  X  *  Xi 

E 

1-15 

Ai  (°K) 

16-30 

Aa  (°K) 

31-45 

Bi  CK/units  of  Xi 

46-60 

Ba  °K/units  of  Xi 

E 

61-75 

*-<P  u  /|i  ),  where  Xi  is  in  feet 

D+2. .£ 

1-15, 16-30, 

Molecular  weights  Mv  for  species 

E 

. .61-75 

1  to  7 

1-15,16-30 

These  values  are:  M].=16.0, 

Ma  =14.0,  M3  =30.0,  Me=32. 0, 

Me  =32.0,  Me =28.0,  M7=30.0 
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V.  DESCRIPTION  OF  OUTPUTS 

The  output  of  the  program  consists  of  a  title  page  con¬ 
taining  program  title,  names  of  originator  and  programmer, 
a  title  statement  describing  the  type  of  computer  run,  date, 
etc.,  and  then  17  lines  listing  the  numerical  values  of  all 
input  data. 

For  each  step  in  the  X  direction  (CSI),  or  horizontal 
coordinate,  results  are  printed  in  a  three  page  format  listing 
the  following  information: 

Page  1  -  The  value  of  "CSI"  followed  by  a  seven 
column  table.  Each  row  of  the  table  represents  data  for  a 
value  of  "PSI,"  or  vertical  coordinate.  The  columns  are, 
from  left  to  right,  (PSI),  stagnation  enthalpy  ratio,  G, 
temperature  (TEMP),  density  (RHO) ,  molecular  weight  of  the 
mixture,  (M) ,  electron  mass  fraction  (ELEC.  CON),  and  elec¬ 
tron  density  in  particles  per  cc. 

Page  2  -  An  eight-column  table  where  the  first 
column  contains  each  value  of  PSI,  while  the  remaining  col¬ 
umns  are  the  mass  fractions  of  each  specie.  The  species  are, 
from  left  to  right,  0,  N,  '  Oi ,  Oa ,  Na ,  and  N0+. 

Page  3  -  A  five-column  table  where  the  first  column 
contains  the  vertical  "ZETA"  variable  corresponding  to  each 
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"PSI"  value.  The  remaining  columns  are,  from  left  to  right, 
incompressible  viscosity  (I-VIS) ,  compressible  viscosity 
(C-VIS) ,  velocity  (U-BAR) ,  and  the  physical  vertical  coordi¬ 
nate  associated  with  PSI  and  ZETA,  the  (Y  coordinate)  in  feet. 

Following  the  data  table  on  page  2  are  printed  two 
integers,  LE  and  LS.  They  indicate  the  number  of  PSI  values 
used  in  the  energy  and  species  solutions,  respectively. 

Following  the  data  table  on  page  3,  except  for  CSI  =  0, 
are  printed  the  values  of  X-coordinate,  PHI,  ZETA  DELTA, 

DDCHI  LOG  SIGMA,  SIGMA  OVER  MU-BAR,  heat  transfer  Q-Dr^  i  ■ 

BTU  per  square  foot- sec,  and  CF . 

Of  the  preceding  quantities,  stagnation  enthalpy, 
density,  temperature,  velocity,  and  viscosity  are  normalized 
with  respect  to  the  input  edge  conditions. 

Examples  of  the  output  described  herein,  appear  in 


Appendix  2. 
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VI .  OPERATING  PROCEDURE 

The  program  was  written  for  the  IBM  709/90/94  digital 
computers  and  uses  the  IBM  FORTRAN  II  monitor  system. 

The  FORTRAN  II  monitor  system  has  standard  tape  desig¬ 
nations,  which  are: 

A2  -  Standard  input  tape 
A3  -  Standard  BCD  output  tape 
Al  -  Systems  tape 

A5  -  Binary  tape  for  restart  procedure. 

An  IOU  subroutine  i»  included  in  the  object  deck  to 
ensure  compatibility  with  the  logical  assignment  of  tapes. 

"Checkpoint"  Procedure 

If  a  restart  option  is  to  be  implemented  a  tape  must  be 
mounted  on  logical  unit  A5.  Depressing  sense  switch  6  at 
any  time  during  the  course  of  a  run  will  dump  the  contents  of 
core  memory  onto  tape  A5  and  then  terminate  the  run.  Tape 
A5  is  dismounted  and  saved  for  future  use.  Tape  A3  may  then 
be  listed. 

To  restart  at  a  future  time,  the  binary  tape  that  was  saved, 
is  again  mounted  on  logical  unit  A5  and  a  small  binary  object 
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deck  labelled  "RESTART,"  is  used  as  the  program  deck.  The 
program  will  read  the  contents  of  tape  A5  anc  processing  will 
commence  from  the  point  where  it  was  formerly  dumped.  Pro¬ 
cessing  will  continue  until  Sense  Switch  6  is  again  depressed 
for  a  second  dump  onto  tape  A5  for  a  future  second  restart, etc. 

To  protect  the  original  information  on  tape  A5,  a  second 
tape  may  be  mounted  on  a  unit  to  be  designated  as  A5  after 
the  original  tape  has  been  read  by  the  7094.  The  unit  with 
the  original  tape  should  be  dialed  off  and  the  tape  dismounted. 
Core  memory  will  be  dumped  on  the  second  tape  for  a  future 
restart . 

The  program  is  normally  terminated  by  specifying  a  value 
of  CSI  FINAL  on  input  card  3.  When  the  program  has  calculated 
the  data  for  the  first  value  of  CSI  which  is  greater  than  CSI 
FINAL,  the  program  will  automatically  process  additional  sets 
of  input  data,  or  in  the  absence  of  such  cards,  will  terminate. 
A  maximum  time  limit  should  be  specified  in  the  instructions 
to  the  operator  in  this  case,  in  the  event  of  a  failure  of  the 
program  to  achieve  a  value  of  CSI  FINAL. 

There  are  several  other  program  stops,  ca;ised  by  numer¬ 
ical  errors,  wherein  the  program  will  print  a  code  number,  and 
in  some  cases  an  alphabetical  statement  describing  the  error. 
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A  list  of  these  error  stops  is  given  in  Appendix  1.  The 
program  will  then  either  process  the  next  data  case,  or 

terminate  just  as  in  the  case  of  a  normal  stop  at  CS1  FINAL. 

Several  options  for  methods  of  numerical  calculation  of 
the  program  can  be  specified  on  input  card  8  and  are  des¬ 
cribed  in  Section  IV.  There  is  one  opcion  controlled  by 
Sense  Switch  1  as  follows: 

SENSE  SWITCH  1 

UP  -  Compressible  viscosities  are  not 

computed  and  printed.  In  this  case, 
either  the  number  0  or  the  values  oi. 
incompressible  viscosity  are  printed, 

the  latter  for  values  of  PSI  greater 

than  PSI  DELTA. 

DOWN  -  Compressible  viscosities  are 
computed  and  printed. 

Sense  Switch  1  instructions  need  be  given  to  the  machine 
operator  only  if  the  Sense  Switch  1  DOWN  option  is  desired. 

See  Eqs.  (15)  and  (16)  for  the  equations  defining 
incompressible  and  compressible  viscosities. 
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NOMENCLATURE 


G 

H 

h. 

T 

r*- 

P 

o 

q 

s 

G 

T 

l 

«k 

X 

z 

Greek  Symbols 
6 

c 


T),  S 


<T  (x)  ,  1?(x)  ,  £(x) 
0 


stagnation  enthalpy  ratio  =  H/H^ 

stagnation  enthalpy  -■  h  +  us/'2 

static  enthalpy  of  species?  h  =  S  h,  Y. 

k  *  k 

effective  Prandtl  number 


hea*  transfer  per  unit  time  per  unit  area 


effective  Schmidt  number 


static  temperature 


mass  averaged  velocity  in  axial  direction 
molecular  weight  of  species  k 


axial  coordinate 


normal  coordinate 


boundary  layer  thickness 

transformed  variable  defined  by  Eq.  (35) 
of  Ref.  1 

transformed  variables  defined  by  Eq.  (53) 
of  Ref.  1 

stretching  function  (see  Eq.  (15)  of  Ref.  1) 
momentum  thickness 
laminar  viscosity  coefficient 
kinematic  viscosity  coefficient 


63 


X 

P 

pe 

y 

<p 

rs* 

St 

Subscripts 

e 


transformed  variable  (see  Eq.  (47))  of  Ref.  1 
mass  density 
eddy  viscosity 
shear  stress 

u /u 

e  T 

stream  function  defined  by  Eq.  (57)  of  Ref.  1 

free  stream 


(") 


incompressible  state 
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LIST  OF  ERROR  STOPS 


PRINTED 

NUMBER 


DESCRIPTION  OF  ERROR 


An  element  in  Column  2  of  species  or  energy- 
difference  equation  matrix  is  equal  to  zero 

The  number  of  ip  values  given  to  Subroutine  HERB 
is  greater  than  149 

A  value  of  ip  greater  than  PSI  DELTA  plus  1/2 
DELTA  PSI  has  been  computed  by  Subroutine  HERB 

Subroutine  HERB  has  taken  more  than  15  itera¬ 
tions  to  compute  a  value  of  ZETA 

No  value  of  'ZETA"  is  greater  than  4^0 

Subroutine  HERB  has  computed  a  value  of  PSI 
DELTA  less  than  1/2  DELTA  PSI 
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